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Abstract 

We study nonlinear dynamics of inhomogeneous DNA double helical chain under dy- 
namic plane-base rotator model by considering angular rotation of bases in a plane 
normal to the helical axis. The DNA dynamics in this case is found to be governed 
by a perturbed sine-Gordon equation while taking into account the interstrand hy- 
drogen bonding energy between bases and the intrastrand inhomogeneous stacking 
energy and by making an analogy with the Heisenberg model of the Hamiltonian 
of an inhomogeneous anisotropic spin ladder with ferromagnetic legs and antiferro- 
magnetic rung coupling. In the homogeneous limit the dynamics is governed by the 
kink-antikink soliton of the sine-Gordon equation which represents the formation 
of open state configuration in DNA double helix. The effect of inhomogeneity in 
stacking energy in the form of localized and periodic variations on the formation of 
open states in DNA is studied under perturbation. The perturbed soliton is obtained 
using a multiple scale soliton perturbation theory by solving the associated linear 
eigen value problem and by constructing the complete set of eigen functions. The 
inhomogeneity in stacking energy is found to modulate the width and speed of the 
soliton depending on the nature of inhomogeneity. Also it introduces fluctuations in 
the form of train of pulses or periodic oscillations in the open state configuration. 
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1 Introduction 



A number of theoretical models have been proposed to describe nonlinear 
molecular excitations in DNA double helix which plays an important role in 
the conservation and transformation of genetic information in biological sys- 
tems [1]. These theoretical models are based on longitudinal and transverse 
motions, as well as bending, stretching and rotations [2,3]. Among the differ- 
ent motions, the rotational motion of bases in DNA is found to contribute 
more towards the opening of base pairs. The first contribution towards non- 
linear dynamics of DNA was made by Englander and his co-workers [4] who 
studied the dynamics of DNA open states by taking into account only the 
rotational motion of nitrogenous bases, which made the main contribution to- 
wards the formation of open states. Yomosa [5,6] developing this idea further 
proposed a dynamic plane base rotator model which is a generalized version 
of the Frenkel-Kontrova [7] model that was later improved by Takeno and 
Homma [8,9] in which attention was paid to the degree of freedom, character- 
izing base rotations in the plane perpendicular to the helical axis around the 
backbone structure. In the above, the DNA dynamics was governed by the 
completely integrable sine-Gordon model admitting kink-type solitons. Then 
Peyrard and Bishop [10] studied the process of denaturation in which only the 
transverse motion of bases along the hydrogen bond was taken into account. 
There was one more model studied by Christiansen and his colleagues (see 
for e.g. [11]) using Toda lattice model in which two types of internal motions 
namely, transverse motion along the hydrogen bond direction and longitudi- 
nal motion along the backbone direction were found to contribute to DNA 
denaturation process in terms of travelling solitary waves and standing waves. 
These localized nonlinear excitations further explain conformation transition 
[12,13,14], long range interaction of kink solitons in the double chain [15,16], 
regulation of transcription [15,17], denaturation [10] and charge transport in 
terms of polarons and bubbles [18]. Some of them have been successfully used 
for interpreting experimental data related to microwave absorption [19,20]. 
Further developments, in this approach for several years was limited to small 
improvements of the models involving only numerical methods of simulation of 
the internal dynamics of DNA [21,22,23]. Also, bubbles [24] discrete breathers 
[25,26,27] and non-breathing compacton-like modes were obtained by solv- 
ing the DNA-lattice model [28]. Eventhough, the rotation of bases in DNA 
is mainly due to thermal forces the thermal fluctuation in DNA dynamics 
has been introduced through random forces in the recent past. For instance, 
Yakushevich et al [23] has shown that topological solitons of the DNA chains 
are stable with respect to thermal oscillations. Since random thermal forces 
introduce only small fluctuations, it is not included in the present study. Thus, 
the study of nonlinear excitations in DNA molecular chain has become an im- 
portant task since it is related to its major functions. 
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In all the above studies, DNA double helix with homogeneous stiff strands has 
been considered for the analysis. However, in nature, the presence of different 
sites along the strands such as promotor, coding, terminator etc. each of which 
having a very specific sequence of bases and particular functions makes the 
strands site-dependent or inhomogeneous(soft). Also, defects caused due to 
external molecules in the sequence and the presence of abasic site-like nonpolar 
mimic of thymine lead to inhomogeneity [29,30]. When included, the DNA 
dynamics is governed by an inhomogeneous perturbed sine-Gordon equation 
and thus the problem boils down to solving the same and finding perturbed 
solitons. Also, in a different context, in the recent times, the study of wave 
propagation, especially solitons through inhomogeneous or disordered media 
assumed lot of interest [31]. For instance, kink-impurity interaction and its 
scattering in the sine-Gordon model was studied in detail by Zhang Fei et al 
[32] . With this in mind, in the present paper we study the nonlinear dynamics 
of DNA double helix with inhomogeneous strands by considering a plane base 
rotator model along the lines of Yomosa. The paper is organized as follows. 
In section 2 we introduce a dynamic plane base rotator model for the site- 
dependent DNA double helix and derive the nonlinear dynamical equation. 
In section 3, we study the effect of inhomogeneity in stacking energy on the 
open state of DNA in terms of kink-antikink solitons by solving a perturbed 
sine-Gordon equation using a multiple scale soliton perturbation theory. The 
results are concluded in section 4. Detailed evaluation of few integrals using 
residue theorem is given in Appendix. 



2 Plane-base rotator Model and Dynamical equation 

In Fig. la we have presented a schematic structure of the B-form DNA double 
helix. Here S and S' represent the two complementary strands in the DNA 
double helix. Each arrow in the figure represents the direction of the base at- 
tached to the strand and the dots between arrows represent the net hydrogen 
bonding effect between the complement bases. While a horizontal projection 
of the n th base pair in the XY-plane is represented in Fig. lb, in Fig.lc, we 
have given the projection of the same in the XZ-plane. The Z-axis is chosen 
along the helical axis of the DNA. In Fig. lb, Q n and Q' n denote the tip of 
the n th bases belonging to the complementary strands S and S'. P n and P' n 
represent the points where the bases in the n th base pair are attached to the 
strands S and S' respectively. Let (9 n , 4> n ) and (9' n , 4>' n ) represent the angles of 
rotation of the bases in the n th base pair around the points P n and P' n in the 
XZ and XY-planes respectively. 

The conformation and stability of DNA double helix is mainly determined 
by the stacking energy between the intrastrand adjacent bases, the hydrogen 
bonding energy between the interstrand complementary bases and other ener- 



3 



(a) 



(c) 



Fig. 1. (a) A schematic structure B-form DNA double helix, (b) A horizontal pro- 
jection of the n th base pair in the XY-plane. (c) A projection of the n th base pair 
in the XZ-plane. 

gies. From a heuristic argument it was assumed that the interstrand base-base 
interaction or hydrogen bonding energy of the given base pair depends on the 
distance between them. Thus, from Fig. lb we can write down the square of 
the distance between the edges of the arrows (Q n Q' n ) 2 as 

(QnQ' n ) 2 = 2 + 4r 2 + (z n - z' n f + 2(z n - z' n ) (cos0 n - cosO 
— 4r [sin 9 n cos n + sin 9' n cos 4>' n ] + 2 [sin 9 n sin 9' n 
x (cos 4> n cos 4>' n + sin 4> n sin (f>' n ) — cos 9 n cos 9' n ] , (1) 

where V is the radius of the circle depicted in Fig. lb. The base-base inter- 
action energy can be understood in a more clear and transparent way by 
introducing quasi-spin operators S n = (S*,Sv,S*) and S„ = (S* , , S„) in 
the form 



S* = sin 9 n cos <f> n , Sl = sin 9 n sin <f> n , S^ = cos 9 n , (2a) 

^ = sin^cos^, ^ = sin^sin0;, S' n z = cos9' n , (2b) 

for the n th bases in the S th and S' th strands respectively. In view of this, 
Eq. (1) can be written in terms of the given spin operators as follows. 
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Fig. 2. A schematic representation of DNA as an anisotropic coupled spin chain 
model or spin ladder. 



(Q n Q> n ) 2 = 2 + 4r 2 + 2 + S*5j - S z n S' n z ] - 4r [s x n + S' n x ] . (3) 

While writing Eq.(3) we have neglected the longitudinal compression along the 
direction of the helical axis thereby choosing z n = z' n . It may be noted that the 
form of {QnQ'n) 2 given in Eq.(3) is the same as the Hamiltonian for a general- 
ized form of the Heisenberg spin model. Therefore, the intrastrand base-base 
interaction in DNA can be written using the same consideration. It is known 
that stacking or base-base interaction is a dominant force that stabilizes the 
DNA double helix. It is much stronger than the hydrogen bonding force and 
in fact the stacking between adjacent bases often contributes more than half 
of the total free energy of the base pairs [33]. Several non-covalent forces in- 
cluding dipole-dipole interaction, van der Waals force etc stabilize stacking in 
DNA. It is also reasonable to think that if such a quasi-spin model can be used 
in this problem, the double strand DNA and the rung-like base pairs can be 
conceived as an anisotropic coupled spin chain model or spin ladder. This is 
schematically presented in Fig. 2. In the figure, S and S' which previously rep- 
resented two strands of the DNA here correspond to two ferromagnetic lattices 
of a spin ladder with antiferromagnetic coupling among the rungs. In the case 
of spin chains each arrow represents the magnetic moment corresponding to 
a group of atoms at that lattice point. Due to antiferromagnetic type of rung 
coupling the arrows in the lattice S and S' are marked anti-parallel to each 
other. Here Z-direction (i.e) the direction of the helical axis is chosen as the 
easy axis of magnetization in the spin chain. Further, the spin-spin exchange 
interaction is restricted, to the nearest neighbours, (i.e) the n th spin is coupled 
to the spins at the (n + l) th and (n — l) th sites. 



With this consideration we use the following Heisenberg model of the Hamil- 
tonian for an anisotropic coupled spin chain model or spin ladder with site- 
dependent or inhomogeneous ferromagnetic-type exchange interaction between 
nearest neighbouring spins in the same lattice (equivalent to coupling among 
bases in the same strand i.e. intrastrand interaction) and antiferromagnetic 
rung-coupling (between bases belonging to the complementary strands i.e. in- 
terstrand interaction) . 
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n 

+A(S z n f + A'(S':) 2 } ■ (4) 

In the above Hamiltonian J and j' correspond to the intrastrand interac- 
tion constant or the stacking energy between the n th base and its nearest 
neighbours in the plane normal to the helical axis in the strands S and S' 
respectively. When K and K' are not equal to J and J' respectively, they in- 
troduce anisotropy in the intrastrand interaction, /j, and rj represent a measure 
of interstrand interaction or hydrogen bonding energy between the bases of 
similar sites in both the strands along the direction of the helical axis and in 
a plane normal to it respectively. Here we have assumed that there exist on 
an average almost uniform interactions A and A that assume positive values 
which are the uniaxial anisotropy coefficients leading to rotation of the bases 
in a plane normal to the helical axis. The quantities f n and f' n in Hamilto- 
nian (4) indicate that the intrastrand stacking energy between bases in the 
S th and S' th strand varies in a specified site dependent fashion which leads 
to inhomogeneity in the DNA double helical chain. In general f n and f n may 
take different values for different base sequences. The inhomogeneity in DNA 
double helix may arise due to any one of the reasons mentioned in the previ- 
ous section. The effect of inhomogeneity in DNA nonlinear dynamics has been 
studied in the past by several authors in various contexts such as variational 
mass density of bases [34] and defect [25] . 



To proceed further we rewrite Hamiltonian (4) in terms of the variables (9 n , (f> n ) 
and (9' n ,(j)' n ) using Eqs.(2) and obtain 



H = J2[~Jfn sin 9 n sin 9 n+1 cos(0 n+ i - n ) - Kf n cos 9 n cos 9 n+1 

n 

-J'f n sin 9> n sin 9' n+1 co^' n+l - <t>' n ) - K% cos 9> n cos 9' n+1 
+i] sin 9 n sin 9' n cos(0„ — 4>' n ) + fi cos 9 n cos 9' n + A cos 9 2 n 

+A'cosC 2 ] • (5) 

The quasi-spin model thus introduced implies that the dynamics of bases in 
DNA can be described by the following equations of motion. 



o J—^IL s ^L^L e> J_^L x> _^_^L (6) 

n s\n9 n d<Pn sin 9 n d9 n J n sin 9< n d^ n ' <?n sin 9' n d9' n 1 ' 

In Eq.(6) the overdot represents time derivative. When the anisotropy energies 
A and A' are much larger than the other interactions, (i.e) when A, A' » 
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J, J', K, K', 77, fj,, then the equations of motion (6) on substituting the Hamil- 
tonian (5) become 



0„ = 2Acos6 n , $ n = 2A'cos9' n . 



(7) 



The other two equations in (6) satisfy identically. Using Eq. (7) in the Hamil- 
tonian (5) we obtain 



J f n sin 9 n sin 9 n+1 cos(0 n+ i - n ) - Kf n cos 9 n 



x cos 9 n+1 - J'f n sin 9' n sin 9' n+1 cos(^ +1 - <t>' n ) - K' f n cos 9> n cos 
+r? sin 9 n sin 6^ cos(0 n - (f)' n ) + [i cos 9 n cos #'J , (8) 

where / = ^ an d ^' — The above Hamiltonian can be rewritten in the 
limits of plane-base rotator {9 n = 9' n = n/2) and absolute minima of potential 

as 



H 



E 

n 

+ J'fn 



I ■ I' •; 



1 - C0S(<^4 +1 - (f>' n ) -7] [1 - COS(0„ - <f>' n )]] . 



(9) 



It may be noted that the above limit corresponds to the XY-spin model of 
two coupled inhomogeneous ferromagnetic spin system. In Hamiltonian (9) 
the first two terms represent the kinetic energies of the rotational motion of 
the n th nucleotide bases accompanied by the potential energy associated with 
the n th nucleotide sugar and phosphate and its complementary unit around 
the axes at P n and P' n (see Fig. lb). / and /' are the moments of inertia of 
the nucleotides around the axes at P n and P' n respectively. It may be further 
noted that in the new Hamiltonian the term proportional to fi vanishes. Now, 
using Hamiltonian (9), the Hamilton's equations of motion can be immediately 
written as 



Hn = J [fn sin(0„ + i 
I'4>' n = J' [ft sin 



71+1 



- <t>n) ~ fn-i sin(0 n - 0„_i)] + i]sm((j) n - <j>' n ), (10a) 
- O - fn-i sm(0^ - + ^sin(0; - n ).(lOb) 



Eqs.(lOa) and (10b) describe the dynamics of DNA in a plane-base rotator 
model at the discrete level while considering the dominant angular rotation 
of the bases in a plane normal to the helical axis and ignoring all other small 
motions of the bases. 



It is expected that in the B-form of DNA double helix the difference in the 
angular rotation of bases with respect to neighbouring bases along the two 
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strands is small [8,9]. Hence we assume that sin(0 n ±i — ra ) ps (0 ra ±i — <Pn) and 
sin(0^ ±1 - 4>' n ) w {4>' n±1 - 4>' n ) in Eqs.(10a, b). Also, as the length of the DNA 
chain is very large involving several thousands of base pairs compared to the 
distance between the neighbouring bases along the strands we make a contin- 
uum approximation as done by several authors in the past [2,3,4,5,6,7,8,9,10] 
which is valid in the long wavelength, low temperature limit by introducing 
two fields of rotational angles, <p n {t) — > <f>(z,t), (f>' n (t) — > (f>'(z,t) and two inho- 
mogeneous stacking fields, /„ — > f(z) and f' n — > f'(z) along with the following 
expansions. 



d(f) a 2 d 2 (f) df a 2 d 2 f 

n±1 = <f){z, t) ± a— + -^-2 ± f n±1 = f(z) ± a — + --^± .(41) 

where 'a' is the lattice parameter along both the strands S and S'. In a similar 
way we write down expansions for (f>' n ±i and /^±i- As the inhomogeneity here 
is site-dependent and associated with the bases themselves, we have chosen 
same lattice parameter 'a' for both n ±i and f n ±i- Under this continuum 
approximation the equations of motion (10a, b) upto 0(a 2 ) is written as 



I(j) tt = Ja 2 [f(z)(f) zz + f z (f) z ] +r]sin((f) - <fr'), (12a) 
JV4 = J' a 2 [f{z)ct>' zz + />'J + 77sin(0' - 0). (12b) 

In Eqs. (12) the suffices t and z represent partial derivatives with respect to 
time t and the spatial variable z respectively. 



In a DNA chain, the two strands are expected to exhibit similar type of macro- 
scopic physical behaviour and hence we assume that 1 = 1', J = J' and / = /'. 
In view of this, Eqs. (12) after rescaling the time variable as i = yf^f-t and 
choosing rj = — for future convenience can be written as 

<t>U = [f( z )(Pzz + fz<Pz] - ^ sin(0 - 0'), (13a) 

1 
2 



+ fM - o sin(0' - 0). (13b) 



It is more convenient to describe the transverse motion of bases in DNA 
strands in terms of the centre of mass co-ordinates. For this, we rewrite 
Eqs. (13a) and (13b) by subtracting and adding them respectively. 

(0 - <f>'h = f(z)(<l> - <P')zz + M - 0% - sin(0 - 0'), (14a) 
(0 + 0')« = f(z)(<P + <P')zz + f z {<P + 0%- (14b) 
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In a different context while studying the magnetoelastic effect induced by 
interaction between two ferromagnetically coupled XY-spin chains in the static 
limit Dandoloff and Saxena [35] obtained similar equations. To commence the 
open state configuration in DNA, the two complementary bases are expected 
to rotate in opposite directions so that = —(f)', and in this case Eq. (14b) 
satisfies identically and Eq. (14a) becomes 



where \1/ = 20. Assuming small inhomogeneity along the strands by choosing 
f(z) = 1 + eg(z) where e is a small parameter and g(z) a measure of the 
inhomogeneity, Eq. (15) can be written as 



Eq. (16) describes the dynamics of bases under rotation in a plane-base rotator 
model of an inhomogeneous DNA double helical chain. When e = 0, Eq. (16) 
reduces to the completely integrable sine-Gordon equation which admits kink 
and antikink-type of soliton solutions and hence we call Eq.(16) as a perturbed 
sine-Gordon equation. The integrable sine-Gordon equation (e = case) was 
originally solved for N-soliton solutions using the most celebrated Inverse Scat- 
tering Transform (1ST) method by Ablowitz and his co-workers [36]. The kink 
and antikink soliton solutions of the integrable sine-Gordon equation (Eq.(16) 
when e = 0) are depicted in Figs. 3a and 3b. The kink-antikink solitons of the 
sine-Gordon equation describe an open state configuration in the DNA dou- 
ble helix. The formation of open state configuration in terms of kink-antikink 
pair in DNA double helical chain is schematically represented in Fig. 3c. In this 
figure the base pairs are found to open locally in the form of kink-antikink 
structure in each strand and propagate along the direction of the helical axis. 



3 Effect of stacking energy Inhomogeneity on the Open State 

3.1 A Perturbation approach 

When inhomogeneity in stacking is present (i.e) when terms proportional to e 
are present in Eq.(16), the inhomogeneity is expected to perturb the kink and 
antikink solitons corresponding to the open state of DNA. One of the most 
powerful techniques in dealing with perturbed soliton equations is the soliton 
perturbation theory which is based on the 1ST method [37,38]. However, as 
the method is very sophisticated it is very difficult to use the same in several 



% = /^ + /^ z -sin*, 



(15) 



#tt-** a + sinW = e[0(z)ttJ z . 



(16) 
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Fig. 3. (a) Kink and (b) antikink soliton solutions of the sine-Gordon equation (Eq. 
(16) when e = ). (c) A sketch of the formation of open state configuration in terms 
of kink-antikink solitons in DNA double helix. 

cases. In view of this, a direct method to study the soliton perturbation was 
first introduced by Gorshkov and Ostrovskii [39] and later many authors used 
different types of direct methods to study soliton perturbation (see for e.g. 
refs. [40,41,42,43,44,45]). The characteristic feature of this method is that the 
perturbed nonlinear equation is linearized by expanding its solution about 
the unperturbed solution and the eigen functions for the operator associated 
with the linearized equation are found out. The complete solution is then 
written in terms of these eigen functions [46,47]. In these methods the basic 
fact that the presence of perturbation not only modifies the shape of the 
soliton by a correction of linear dispersion but also undergoes a slow time 
change of the soliton parameters have been acknowledged. In this paper, we 
use one such direct perturbation method which is also dealt in reference [48] 
in a different context to solve the perturbed sine-Gordon equation (16) and 
to understand the effect of inhomogeneity in stacking energy on the open 
state configuration of DNA. The procedure we adapt here is based on the 
derivative expansion method to linearize the perturbed sine-Gordon equation 
in the co-ordinate frame attached to the moving frame. The parameters of the 
kink-antikink soliton are assumed to depend on a slow time scale in order to 
eliminate the secular terms. The linearized equations will be solved using the 
method of separation of variables which ultimately will be related to an eigen 
value problem, the eigen functions of which form the bases of the perturbed 
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solution. The eigen functions contain information about the time dependence 
of the soliton parameters and help to calculate the perturbed soliton. In the 
following we use the above approach to find the perturbed soliton solution of 
Eq.(16). 



3.2 Linearization of the perturbed sine-Gordon equation 



When the perturbation is absent (i.e) when e = in Eq. (16) the unper- 
turbed integrable sine-Gordon equation provides N-soliton solutions and the 
one soliton solution (see also Figs. 3a and 3b is written as 



ty(z,i) = 4arct&nexp[±m(z — vi)], m = —= (17) 

V 1 — v 2 

In Eq.(17) while the upper sign corresponds to kink soliton, the lower sign 
represents the antikink soliton. Here v and vrT x are real parameters that de- 
termine the velocity and width of the soliton respectively. In order to study 
the effect of perturbation the time variable t is transformed into several vari- 
ables as t n = e n t where n=0,l,2,... and e is a very small parameter. In view of 
this, the time derivative in Eq. (16) is replaced by the expansion 



d d d 9 d j . 

^ = ^- + e 7 - + e 2 — + .... 18 
at ot oti ot 2 

Simultaneously \I> is expanded in an asymptotic series as 



^ = ^(°) +e ^( 1 ) + eV 2) + .... (19) 

Using the above expansions for t and ^ in Eq. (16) and equating the coefficients 
of different powers of e, we obtain the following equations. 



e (0) : <l-^i°2 + sm^ (0) =0, (20a) 
■■ *$ - *S + cos = g*V> + g z mf - 2^ , (20b) 

etc. 

The initial conditions for perturbation in the case of the single soliton given 
in Eq. (17) is written as 



^ (0) (^,0) = 4arctanexpm^,^ (ri) (^,0) = * t ( " } (^,0) = 0, n = 0,1, ....(21; 
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The equation obtained at order of e^, (i.e) Eq.(20a) is just the integrable sine- 
Gordon equation for tp(°\ the single soliton solution of which can be written 
from Eq. (17) immediately as 

¥°\z,t ) = 4arctanexpC, C = ±m (z - £), & = v , (22) 

where v is the velocity of soliton in the to-time scale. Due to perturbation, the 
soliton parameters namely m and £ are now treated as functions of the slow 
time variables to,ti,t2, •••• However m is treated independent of to- m view of 
the above, Eq. (20b) becomes 

*$ - *£? + (! - 2sech 2 ()^ = F W(z), (23a) 
where 



F^\z) = 2 [g(z)sechz} z + Avosechz m tl + {m 2 C, tl — zm tl ) tanhz . (23b) 



While writing the above equation we have used the result cos\l/^ ^ = 1 — 
2sech 2 ( obtained from the solution of the unperturbed equation (20a). In 
order to represent everything in a co-ordinate system moving with the soliton, 
we transform z = ^ + vt and t = t + (, so that Eqs. (23) and the initial 
conditions given in Eq. (21a) can be written as 



vJ>$ - 2mv ^ - + (1 - 2sech 2 ()^ = F^((, t ), (24a) 
where 



F«(C,t ) = 2 [g(()sech(] c + 4v sech( [m tl + (m 2 ^ - (m tl ) tanhC](24b) 

and 



¥ 1 \(,0) = 0, tfg>(C,0) = 0. (24c) 

We now introduce one more transformation r = J^- — on the indepen- 

dent variable to eliminate the first term in the left hand side of Eq. (24a). 
Thus on using the above transformation Eq. (24a) becomes 



M>« - + (1 - sech 2 ()¥V = F«(C,t), (25) 

where r) equals exactly the right hand side of Eq. (24b). The solution 

of Eq. (25) is searched by assuming 
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* (1) (C,r) =X(C)T(r), FV((,t)=X c (()H(t). 
Substituting Eq. (26) in Eq. (25) we obtain 



(26) 



Xr 



X« + (2sech 2 ( -l)X = - [T T - H{r)\ . 



(27) 



In Eq. (27) the left hand side is independent of r and the right hand side is 
independent of the variable (. Hence we can equate the left and right hand 
sides of Eq. (27) to a constant say Aq and write 



X cc + (2sech 2 ( - l)X = X X o (28a) 
T T - A T = H(t). (28b) 

Thus, the problem of constructing the perturbed soliton at this moment turns 
out to be solving Eqs. (28a) and (28b) by constructing the eigen functions and 
finding the eigen values. It may be noted that Eq. (28a) is a generalized eigen 
value problem which is not a self- adjoint eigen value problem and differing 
from the normal eigen value problem with X^ in the right hand side instead 
of X and Eq. (28b) is a first order linear inhomogeneous differential equation 
which can be solved using known procedure. 



3.3 Solving the eigen value problem 



Before actually solving the eigen value equation (28a) we first consider it in a 
more general form given by 



L±X = XX, Li = % + 2sech 2 ( - 1, (29) 

where A is the eigen value. In order to find the adjoint eigen function to X, 
we consider another eigen value problem 



L 2 X = AX, (30) 

where the operator L 2 is still to be determined. The eigen value equations (29) 
and (30) can be combined to give 



L 2 LiX = X 2 X, L X L 2 X = X 2 X. (31) 

Since we already know the form of L\ as given in Eq.(29), if we choose the 
operator L 2 as L 2 = + 6sech 2 ( — 1, it can be verified that LiL 2 is the 
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adjoint of L 2 Li. Thus X and X are expected to be adjoint eigen functions. 



Now for solving the eigen value equations (29) and (30), we choose the eigen 
functions X and X to be in the form 



X((,k)=p((,ky k ^ X((,k) = q((,ky h ^ (32) 

where p((, k) and q((, k) are assumed to have the asymptotic behaviour p((, k) - 
a constant and q((, k) — > a constant as ( — > ±oo and fc is the propagation 
constant. On using these asymptotic forms for p((, k) and q((, k) in Eq. (32) 
and then substituting the resultant X((, k) and X((, k) in Eqs. (29) and (30) 
we obtain the eigen value as 



A = -(A: 2 + 1). (33) 

On substituting the exact forms of X and X from Eq. (32) in Eqs. (29) and 
(30) we get the following set of ordinary differential equations for p((, k) and 

?(£*)■ 



(Li - k 2 )p + 2ikp^ + (1 + k 2 )q = 0, (34a) 
(L 2 - k 2 )q + 2ikq c + (1 + k 2 )p = 0. (34b) 

For solving the above equations we expand p((, k) and q((, k) in the following 
series [48]. 



IN sinh( 1 sinh( 1 

P(C,«)=Po +Pi TT +P2 T27 +P3 T3T +P4 rj- + 35a 

coshC cosh C cosh d C cosh ( 

sinhC 1 sinhC 1 N 

9(C ' fc) = 90 + 9l oo^c + g2 c^si7c + g3 oo^i?c + g4 oo^c-' (35b) 

where the coefficients pj and g^, j=0,l,2,... are functions of k which are to 
be determined. We substitute the series expansions given in Eqs. (35a) and 
(35b) in Eqs. (34) and collect the coefficients of 1, eos ^ and obtain 
the following algebraic equations. 



Po = Qo, Pi = Qi, (36a) 

2p + 2ik Pl + (3 - k 2 )p 2 - 4ikp 3 + (1 + k 2 )q 2 = 0, (36b) 

6q + 2ik qi + (3 - k 2 )q 2 - Aikq 3 + (1 + k 2 )p 2 = 0, (36c) 

-Aikp 2 + (3 - k 2 )p 3 + (1 + k 2 )q 3 = 0, (36d) 

4gi - 4iA;g 2 + (3 - A; 2 )g 3 + (1 + k 2 )p 3 = 0. (36e) 
etc. 
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By assuming pj = qj = for j > 3, and substituting these values in Eqs.(36d) 
and (36e), we obtain 



P2 = 0, q 2 = -ij. (37) 
On substituting the results given in Eq.(37) in Eqs. (36b) and (36c) we get 



Pl = qi = -JT^¥-y q2 = -JT^¥-y (38) 

For our convenience we choose po = go = c(l — k 2 ) and on substituting this in 
the above equations, we obtain the other coefficients as follows. 



Pi = qi = -2ik, p 2 = 0, q 2 = -2c. (39) 

Here 'c' is an arbitrary constant which will be determined. Using the above 
values in Eqs. (35a) and (35b), and then in Eq. (32), we finally obtain 



X(C, k) = c(l - k 2 - 2ik tanh ()e ik<: , (40a) 
X(C, k) = c(l - k 2 - 2ik tanh C - 2sech 2 ()e ik<: . (40b) 

On comparing Eqs. (40a) and (40b) we can write 



*«,*) = (41) 

Using Eq.(41) in the right hand side of Eq.(29) and comparing the resultant 
equation with (28a), we can write down the eigen value Aq as 



A. = £±^. (42) 



Now to determine the constant 'c' we use the orthonormality relation between 
X(C, k) and X((, k) given by 

/oo _ 
X((,k)X*((,k')d( = 6(k-k>). (43) 
-oo 

On substituting the eigen functions X((,k) and X((,k) given in Eqs. (40a) 
and (40b) respectively in Eq.(43) and after evaluating the integral we obtain 
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27TC 2 



(1 - fc 2 )(l - k' 2 ) + 4kk'} 6(k - k') = 5(k - k'), (44) 



which gives c 2 = + k 2 ) 2 . The correct form of the eigen functions X(£, k) 
and X((, k) is written down after using the above value of 'c' in Eqs. (40a) 
and (40b). 



(1 - k 2 - 2i/ctanhC - 2sech 2 () ikt 
X(C, fc) = ^= —e lk \ (45b) 



It can be verified that the operator L 2 also has the following two discrete eigen 
functions. 



X (C) = (1 -CtanhC)sec/iC, ^i(C) = sech(tanh (. (46) 

It may be further noted that the eigen function X± corresponds to the discrete 
eigen value A = 0. That is 



L 2 X 1 (C) = 0. (47) 
3.4 Complete set of orthonormal basis 



Having found the eigen functions X((, k) and X((, k), we now check the com- 
pleteness of them by writing 

/oo _ 
X((, k)X*(C, k)dk + /(c, CO = S(( - ('), (48) 
-00 

where /(£, (') is an arbitrary function to be determined. First we evaluate the 
integral in the left hand side of Eq. (48) after substituting the values of the 
eigen functions X((,k) and X((,k) given in Eqs. (45a) and (45b). Thus we 
have the following integrals to evaluate. 

f 00 x(c, k)x*((>, k)dk = 5(( - - - f r T^m em ~ C) 

J-00 7T |_J-oo (1 + At) 

x{2 - sech 2 C - iA;(tanhC - tanhC')} 

x(l -ifctanhC) (ifc + tanhC')]- (49) 
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The integrals in the right hand side of Eq.(49) are evaluated with the aid 
of the residue theorem. It may be noted that the integrands in these two 
integrals as functions of the complex variable k are analytic everywhere in 
the complex A;-plane except at the two poles k — ±i of first and second or- 
der respectively. Let R\ and R 2 be the residues corresponding to the func- 
tions ( T ^e^-^{2-iA;(tanhC-tanhC / ) - sech 2 ('} and -^^e^'^il - 
ik tanh ()(ik + tanh (') tanh (' in Eq. (49) at the pole k = i of first and sec- 
ond order respectively. On calculating the residues Ri and R 2 using standard 
procedure we obtain 



i?i = — [2 + tanh C - tanh (' - sech 2 (] e (c_c,) , (50a) 

R2 = ^ [(1 - C + CO (tanh C - tanh (') - (( - (') 

x(l - tanh C tanh C')] e (c ~°. (50b) 
On summing up the above results Eq. (49) becomes 

^ X((,k)X*((',k)dk = 5((-C) - [sech(sech('(l -C'tanhC') 
J— 00 

+(sech(sechC tanh ('} , (51) 

which can be identified as 

/oo _ r _ , 

x((, k)x*(C, k)dk = 5(( - CO - x (C)x (CO + *i(C)*i(0 , (52) 
-00 

upon using the values of X (Q and X 1 (() as given in Eq.(46). By substituting 
Eq.(52) in Eq.(48), we can evaluate the following value of f((, (') in terms of 
the discrete orthogonal states. 



/(CO = E .Y,(0-YA)- (53) 
3=0,1 

Thus, by comparing Eqs.(51) and (52), we can write two additional orthogonal 
discrete states given by 



X (() = sech(, X 1 (()=(sech(. (54) 

It may be checked that the following relations exist between the discrete states 
X (C),Xo(C),*i(C) and X X (C). 



LiXo = 0, LiXiiC) = -2Xi (C), L 2 X (C) = 2X (C), (55) 
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in addition to the relation given in Eq.(47). In conclusion, we have a set of two 
complete orthonormal bases {X} and {X} given by {X(£,k),Xo(Q,Xi(Q} 
and {X((, k), X ((), X^Q} respectively. These set of orthonormal basis func- 
tions will be used to construct the perturbed soliton solution. 



3.5 Evaluation o/T(r) 



Having solved Eq.(28a) by finding X(Q, in order to construct V^(Cj r), we 
now find T(r) by solving Eq. (28b). For this first we rewrite Eq. (28b) by 
replacing the function H(t) in the right hand side using Eq. (26) and (41). 



T T - A T 



ikX((,k) 



(56) 



On multiplying and dividing the right hand side of Eq. (56) by X*((, k) and 
on integrating with respect to ( between the limits — oo to +oo we obtain 
due to orthonormality of the functions X and X (see Eq.(43)) the following 
equation. 



T T - \ T — — j ^ F«(C, t)X*((, k)d(, 



(57) 



which can be explicitly written after substituting the values of FW((,t) and 
X*((,k) from Eqs. (24b) and (45a) respectively as 



T - _ A ° T = • /o-Jn ^u2\ r dC ([9((>ech(] c + 2v sech( [m tl 
+ (m 2 £ tl -Cm tl )tanhC]) (1 - k 2 + 2ik tanhC)e _ifcC . 



(58) 



On solving the above equation using standard procedure, we obtain 



T(r, k) = C(k)e x " T - 



I d( ([g(()sech(L + 2v sech( 



X 



iV27r\ k(l + k 2 ) 
m tl + (m 2 £ tl -Cm tl )tanhC]) (1 - k 2 + 2ikta,nh()e~ ik( > , (59) 



where C(k) is a constant which can be found by using the initial condition 
T(r, k) = when r - 



Thus we obtain 
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C(k) 



and hence 



. , , , M , ([g(()sech(] c + 2v sech( [m tl 

zV2n\ok(l + k l ) J-oo v s 

+ (™ 2 £ tl - (m tl ) tanhC]) (1 - k 2 + 2ik tanhC)e Ao(1+U(,) 2- i ^, (60) 



,Ao[r+^C'] 



T(r, fc) = / dC ({g(C)sechC} , + 2v sech(' [m tl 

\ / 2m\ k(l + k ) J-oo v s 

+ (m 2 ^ - C'mtJtanhC']) (1 - k 2 + 2^tanhC')e~ ifeC ' 

l) . (61) 



x e 



The allowed discrete values of T will be determined in the next section while 
constructing the perturbed part of the soliton ^W((,t) using the values of 
X(Q and T(r). 



3. 5 Perturbation of soliton 



As mentioned, we now write down the first order perturbation correction 
^ (1) (C,r) in terms of the basis functions {X} = {X (( , k) , X (() , X^Q} and 
{T} = {T(r, fc), T (r), Ti(r)}. As per Eq.(26), in terms of the basis functions 
the perturbed part of the soliton can be written as 



*W(C,t) = / X(C, k)T(r, k)dk + £ ^(C)T,(r). (62) 
J -°° j=o,i 

However, it should be noted that in the basis {T}, T (r) and Ti(r) are yet to be 
determined. Therefore before evaluating the values of ^^(C, k), we determine 
T (t) and Ti(r). This is carried out by substituting Eq.(62) in Eq.(25) which 
finally becomes 



/oo 
ik [T T (r, fc) - A T(r, fc)] X(C, A;)rfA; + T lr (r)X (C) 
-oo 

-[T 0r (r)-2T 1 (r)]X 1 (C)=F( 1 )(C,r), (63) 

upon using Eqs.(29),(41),(47) and (55). Now multiplying Eq.(63) by X*((, k),X (() 

and -Xi(C) separately and using the orthonormal relations such as Jf^ X(C, /c)dA; = 

i,/^x(c,fc)A,-(c)dC = Ho^CC, ^)x,(C)rfC = o > n o x i (c)x,(c)dC = <W, J = 
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0, 1, we get 



; poo 

T t (t, k) - A T(r, k) = - F«(C, t)X*((, k)d(, (64a) 

%K J —oo 

/oo 
FW(C,r)X (C)dC, (64b) 

/oo 
F( 1 )(C,r)X 1 (C)dC- (64c) 
-oo 

As F^((,t) given in Eq. (24b) does not contain time V explicitly, the right 
hand side of Eqs. (64a-c) also should be independent of time. Then it is easy 
to verify that for the values of F^((, r), X*((, k),T(r, k) and A given respec- 
tively in Eqs. (24b), (45a), (61) and (42), Eq. (64a) satisfies identically. As the 
right hand sides of Eqs. (64b) and (64c) are also independent of time, they 
give rise to secularities and hence the nonsecular conditions can be written as 

/° O F( 1 )(C,r)X (C)rfC = 0, (65a) 

J — oo 

/oo 
F( 1 )(C,r)X 1 (C)rfC = 0. (65b) 
-oo 

Using Eqs. (65a) and (65b) back in Eqs. (64b) and (64c), we obtain Ti(t) = 
and T (t) = C\ which has to be determined. For this, we substitute Ti(t) = 
in Eq. (64c) and integrate with respect to r to obtain 

T (r) = d = {l + v) l°° d( ([g(()sech(l + 2v sech( 

J —oo ^ 

X 



m h + (m £ tl - (m h ) tanh( J ( sech(. (66) 



3. 7 Variation of soliton parameters 



We now estimate the nonsecularity conditions (65a) and (65b) by evaluat- 
ing the integrals after substituting the values of F^((, r), Xq(Q and ATi(C) 
respectively from Eqs. (24b) and (54). The results give 



1 f°° 

m tl = -— [g(()sech(} c sech(d(, (67a) 

Z,Vq J —oo 

61 = ~2^k: JZ {9((>ech(] c (sech(d(. (67b) 



Eq. (67a) describes the time evolution of the inverse of the width of the soliton 
and Eq. (67b) gives the velocity of the soliton. g(Q that appears in the above 
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nonsecularity relation is related to the inhomogeneity in stacking energy of 
DNA. In order to evaluate the integrals in Eqs. (67a) and (67b) explicitly 
we have to substitute specific value for g((). We consider g(() in the form of 
localized and periodic functions. A localized g(() corresponds to the interca- 
lation of a compound between neighbouring base pairs without disturbing the 
base pairs and their sequence in the DNA double helical chain. The periodic 
nature of g(() represents a periodic repetition of similar base pairs along the 
helical chain. We consider the localized form of g(() as (i) g(() = sech( and 
the periodic form of g(() as (ii) g(() = cos (. We substitute the above values of 
g(() one by one in Eqs. (67a) and (67b) and evaluate the integrals in the right 
hand side to understand the time evolution of the width of the soliton and its 
velocity. At this point it is worth mentioning that Dandoloff and Saxena [35] 
realized that in the case of XY-spin chains the model of which identifies with 
our plane-base rotator model, the ansatz g(() = sech( energetically favours 
the deformation of spin chains. 



When we substitute g(() = sech( in Eqs. (67a) and (67b) and on evaluating 
the integrals, we obtain 



m " = °- 6. = taK£- < 68) 

The above equations can be rewritten in terms of the original time variable t 
by using the transformation Jt = + e-J^- or in otherwords m t - = m to + em tl 
and £f = £ to + e £ti- As m is independent of to(m to = 0) and £ to = t> , we can 
write 



m = m , ^i = v = v + 2 , (69) 

where 1/mo is the initial width of the soliton. The first of Eq.(69) says that 
when g(() = sech(, the width (m _1 ) of the soliton remains constant. However 
from the second of Eq.(69), we find that the velocity of the soliton gets a 
correction. As the correction term is a definite positive quantity the velocity 
of the soliton increases in this case. Interestingly, the amount of increment in 
velocity depends on the initial width and initial velocity of the soliton. Wider 
the soliton, greater the increment in velocity due to inhomogeneity. Also slowly 
moving solitons gain more speed. The increase in speed helps to overcome the 
barrier of the local inhomogeneity (which may be due to the presence of abasic 
site or intercalation of a molecule) and the solitons representing the open state 
will propagate easily along the chain without formation of a bound state. In a 
similar study on resonant kink-impurity interaction and kink scattering in the 
sine-Gordon model Zhang Fei et al. [32] observed that if the initial velocity of 
the kink is smaller than a critical velocity it will be either trapped or reflected 
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by the impurity. In fact they have showed that for most of the initial veloci- 
ties, the kink is trapped except in the case of some special initial velocities the 
kink may be totally reflected by the impurity. It was further found that when 
the kink velocity is greater than the critical velocity, it will pass through the 
impurity. It is interesting to note that our results on the velocity of soliton 
found in Eq. (69) is similar to the last case of Zhang Fei et al [32] where the 
soliton will pass through by overcoming the barrier of the local inhomogeneity. 

Next, we substitute the periodic function g(Q = cos( in Eqs.(67a) and (67b) 
and evaluate the integrals to obtain 



From the above the parameters m and £ can be written in terms of the original 
variable t after solving Eq.(70) as 



From Eq. (71) we find that the width of the soliton in this case remains 
constant and the soliton velocity increases similar to the case g(() = sechC,. 
On comparing Eqs. (69) and (71), we observe that the increase in velocity of 
the soliton is less in this case. This is because in this case, the inhomogeneity 
occurs periodically in the entire DNA chain in terms of sequence. 



3. 8 First order perturbed soliton 



Now, we explicitly construct the first order perturbation correction to the one 
soliton for the different cases of g(() by substituting the values of the basis 



functions {X} = {X((, k), X ((), X^)} and {T} = {T(r, k), T (r), T^r)} 



after using the values of F^((, r), m tl and £ tl for the respective g(() values 
in Eq. (62). Thus in the case of g(() = sech(, we substitute the values of 
X((,k),X (() and X^Q from Eqs. (45a) and (54) and T(r, k), T (r) from 
Eqs. (61) and (66) and F«(C,r) from Eq. (24b) and use the values of m tl 
and £ tl from Eqs. (68) in Eq. (62) to obtain 



m tl = 0, &i 




(70) 



m = m , £ { = 



v = v + 




(71) 
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x(l — k 2 + 2ik tanh £) (n — 6sec/iC)sec/i£tanh£ 



x 



(72) 



where a = to m j- 1+v °^ an d 8 = (l±gg)(l±*: ) _ ^ While writing the above, 
we have also used Ti(r) = 0,r = ^[^o ~ m (l + v o)C] an d the eigen value 
Ao = %( ~ 1+ k k - ■ The integrals in the right hand side of Eq.(72) can be evaluated 
using standard residue theorem . The details are given in Appendix-A. The 
final form of \I/^(C, to) after evaluating the integrals becomes 



^(C^-^V^C e-^ + ^tanhCy/^C 



7Y 



6v 2 



[2u(l + u)C + + 4cw - l] sech(. (73) 



While constructing the above perturbed solution, we have used the values of 
m's and £'s and of course the corresponding to) values. Finally, the per- 

turbed one soliton solution that is ^(z, t ) = ty(°\z, t ) + ^^\z, t ) (choosing 
e = 1) is written using Eqs. (22) and (73) as 



ty(z, to) ~ 4arctanexp[±m (-2 — i>oto)] + r^j^ \J sec ^ l V L - m { z ~ v ^o)] 

3(m(z-vt ) 64 .,/ — " — 

xe + 2 H ^= tanh[±m(z — t> t )J y secn[±m{z — vto)\ 

27y 2 



x e =F§m(z-yto) _|_ 



6mn 2 



m(iT - 1) + 2t w sec/i[±m(^ - vt )\. (74) 



In Eq.(74) while the upper sign corresponds to perturbed kink-soliton the 
lower sign represents the perturbed antikink-soliton. The rotation of bases 
denoted by <p(z, to) can be immediately written down by using the relation 
<f) — Tf. In Figs. 4a and 4b, we plot 0(^,t o ) (rotation of bases under pertur- 
bation) for the parametric choice Vq = 0.4. From the figure, we observe that 
there appears fluctuation in the form of a train of pulses closely resembling 
the shape of the inhomogeneity profile in the width of the soliton as time 
progresses. Also, as time passes, the amplitude of the pulses generating this 
fluctuation increases. However, in the asymptotic region of the soliton there 
is no change in the topological character and no fluctuations appear in that 
region. It shows that the localized inhomogeneity in stacking energy in DNA 
in the form of a pulse (g(z) = sechz) does not affect very much the opening 
of bases except fluctuations in the form of a train of pulses in the localized 
region of the kink and antikink-soliton. We have schematically represented 
this in Fig. 4c. 
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Fig. 4. (a) The perturbed kink-soliton and (b) the perturbed antikink-soliton for 
the inhomogeneity g(z) = sechz and vq = 0.4. (c) A sketch of the base pair opening 
in DNA double helix with fluctuation. 

We then repeat the procedure for constructing the perturbed one-soliton so- 
lution in the case of g(() = cos( which is written as 

1 poo rjh poo 

* (1) (C, *))«-/ 77-7^(1 - k 2 - 2ik tanhC)e^ / dC(l - k 2 

71 J-oo 1 + K) J-oo 



- x \ 1 + k 2 1 

TT 2 



+2ik tanh () [sin ( + (cos £ ) tanh (}sech( 

8 

x{e i ^ a e iK - e~ iK ] + (1 + v)sech( 



oo ^2 



x/ c?C C 2 [sin C + (cos C- tanh (}sech 2 (. (75) 



oo 



The details of values of the integrals in the above equation are given in 
Appendix-B. The perturbed kink (upper sign)-antikink (lower sign) one soli- 
ton solution in this case is finally obtained as 



7T 2 

ty(z, to) ~ 4 arc tan exp[±m (z — fo^o)] 



m(v 2 - 1) + 2vt 



16mv 2 

xsech[m(z — vto)]. (76) 
After finding <f>(z, to) from the relation = y we plot it in Figs. 5a and 5b for 
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l ining a Mim m 



Antikink 

(C) 

Fig. 5. The perturbed (a) kink-soliton and (b) antikink-soliton for the inhomogene- 
ity g(z) = cosz with vq = 0.4. (c) A sketch of the open state configuration in DNA 
with small fluctuations. 

the same value of the parameter as before. From the figures we observe that 
periodic oscillations appear in the width of the soliton without any change 
asymptotically. 



4 Conclusions 



In this paper, we studied the nonlinear dynamics of DNA double helix with 
stacking inhomogeneity by considering the dynamic plane-base rotator model. 
The dynamical equation which finally appeared in the form of a perturbed 
sine-Gordon equation was derived from a suitable Hamiltonian in analogy 
with Heisenberg model of an inhomogeneous anisotropic coupled spin chain 
or spin ladder with ferromagnetic legs and antiferromagnetic rung coupling in 
the continuum limit. In the unperturbed limit which is also the homogeneous 
limit, the dynamics is governed by the kink-antikink soliton of the integrable 
sine-Gordon equation which represents the open state configuration of base 
pairs in DNA double helix. Even though DNA double helix is a large molecular 
chain involving very large number of base pairs, base pair opening is limited 
to a very few number of base pairs forming localized coherent structure in the 
form of kink-antikink solitons which is obtained as a balance between disper- 
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sion and nonlinearity traveling with constant speed and amplitude without 
loosing its energy along the helical chain. To understand the effect of stacking 
energy inhomogeneity on the open state configuration of base pairs we carried 
out a multiple scale soliton perturbation analysis. For implementing this we 
linearized the perturbed sine-Gordon equation using multiple-scale expansion 
to obtain linearized equations in the form of eigen value problems. The per- 
turbed kink-antikink soliton solutions were constructed for different forms of 
inhomogeneities by solving the associated eigen value problem. By using the 
complete set of eigen functions thus obtained as the basis functions the per- 
turbed solutions were constructed. 



The perturbation not only modifies the shape of the soliton but also undergoes 
a slow time change of the soliton parameters namely the width and velocity 
for the different forms of the inhomogeneity chosen. We chose the inhomo- 
geneity in the form of localized and periodic functions. The results show that 
when the inhomogeneity is either in the form g(z) = sechz or g(z) = cos z, 
the width of the soliton remains constant. Thus, the number of base pairs 
participating in the opening do not change due to the above pattern of inho- 
mogeneities. However, in this case the speed of the soliton increases with a 
correction that is proportional to the square of the initial width of the soli- 
ton and inversely proportional to its initial velocity. Thus the inhomogeneity 
increases the speed with which the base pairs are opening and closing or wind- 
ing and unwinding. From the perturbed solutions corresponding to different 
inhomogeneities (see Figs. 4,5) we observe that the perturbation due to in- 
homogeneity in stacking energy along the strands introduces fluctuation only 
in the width of the solitons. The nature of the fluctuation varies depend- 
ing on the type of inhomogeneity. In particular, when the inhomogeneity is 
in the form of g(z) = sechz, we find that fluctuation in the form of pulse 
trains resembling the shape of the inhomogeneity is generated in the width 
of the soliton representing open state configuration. It is noted that in the 
long time limit, eventhough neighbouring pulses overlap the train of pulse-like 
fluctuations maintain their character without affecting the soliton. In a simi- 
lar way, when the inhomogeneity is of the form g(z) = cos z, the fluctuation 
appears in the form of periodic oscillations in the width of the soliton. In all 
these cases, asymptotically the kink-antikink soliton shape is preserved. The 
results indicate that inhomogeneity in stacking energy in DNA double helix 
can (i) introduce small fluctuations which may merge asymptotically during 
the process of opening and closing of base pairs (ii) increase or decrease the 
number of base pairs participating in the open state configuration and (iii) 
change the speed with which the open state configuration can travel along the 
double helical chain. Thus in conclusion the inhomogeneity in stacking does 
not affect the general pattern of base pair opening. Even though the size of 
the base pair is big and the solvent effect on pure rotation of base pairs is 
negligible, the soliton, a coherent structure which is formed by involving few 
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base pairs move along the helical chain without dissipation or any other form 
of deformation. Similar conclusion was also arrived in the case of propagation 
of soliton representing base pair opening in a discrete site-dependent DNA 
[21] and propagation of bubble in a heterogeneous DNA chain [24]. As nature 
selects generally inhomogeneous DNA, the functions such as replication and 
transcription can be explained more viably through formation of open states 
through our inhomogeneous model rather than the homogeneous case. This is 
also because, it is known that transcription and replication are sequence de- 
pendent. This is further similar to what was observed in the case of proteins 
where inhomogeneity of the sequence leads to inhomogeneous fluctuations, en- 
hanced by the nonlinear effect [49]. Eventhough the relevance of these effects in 
the DNA to biological processes are not yet clearly, established a recent study 
suggests that thermally induced base pair opening agrees with experimental 
observation on DNA base pair opening detected by potassium permanganate 
foot printing [50]. We will make numerical analysis of the discrete dynam- 
ical equation separately to understand the discreteness effect which will be 
published elsewhere. It is also equally important to analyse the nonlinear dy- 
namics of DNA double helix when the hydrogen bonding energy depends on 
the distance between the bases (inhomogeneity in hydrogen bonds) and the 
study is under progress. Also, the nonlinear dynamics study of open state 
configuration facilitated by enzymes (protein) is under progress. 
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A Evaluation of integrals in Eq.(72) using residue theorem 

In this appendix we evaluate the integrals found in Eq.(72) using standard 
residue theorem. Eq.(72) is of the form 



* (1) (C,*o) 



I rOO 



J—oo 



POO 

/ dC(l - *? 




(A.I) 
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The evaluation of various integrals in Eq. (A.l) can be facilitated by first 
finding the values of the following two integrals [51]. 



/oo 
sechC e lxC d(, (A.2) 
-oo 



tanhCe lxC dC, (A.3) 

where x can fake values —k and (5. The integrand in I\ is found to be analytic 
everywhere except at the pole ( — > i(2n + 1)| , n = 0, 1, 2, .... The residue of 
the function sec/i^ e lxC - can be written as 



7T °° P*^ 

ites(C = i(2n + l)-) = V Zim -r- — — , (A.4) 

V ^ of ^(2n + l) f |(cOShC)' 

which can be simplified to give 



Res(( = i(2n + 1)-) = — — . (A.5) 

vs v J 2 J 2% cosh fx 

Thus the integral value of h in Eq. (A.2) is written as 



/oo -77- 
sech^dC, = — ; — t-= — r. (A.6) 
s s cosh (fx) V ; 

Similarly, using the same procedure we evaluate the integral I2 in which the 
integrand tanh£ e 1 ^ contains poles at ( — > ?(2n + l)f, n = 0,1,2,... and 
obtain 



/OO 777- 

Now using the value of the integral I\ as given in Eq. (A.6), we evaluate the 
following integrals found in Eq. (A.l) by integrating them by parts successively 
to obtain 



JZ sec ^ nhCe ' x<dC= ^£w)' (A ' 8a) 



/. 



tt(1 - x 2 ) 



sech( tanh 2 C e ixC d( = v - * \ . (A.8b) 

2 cosh (fx) 1 ' 
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Similarly, using Eq. (A. 7) we evaluate the following integrals found in Eq. (A.l) 
by making successive integrations by parts and obtain 



/. 



rfCtanhC e^d( = - (A.9a) 
oo 2smh(|x) 

sech 2 ( tanh 2 C e ix ^d( = X \ - X ( . (A.9b) 

6sinh(f X ) 1 ; 



Here also x can take values —k and (3. Substituting the values of the integrals 
found in Eqs. (A. 8a), (A. 8b) (A.9a) and (A.9b) in Eq. (A.l), we obtain 



* (1) (C,to) = - 



12 



j — i 



7T 



dk- 



[1-k 2 -2ifctanhC) 
k 2 (l + k 2 ) 2 



- v 2 )k(l + k 2 ) 



xsech'-l3 - {(2 -v)(l + v) 2 (l + k 2 ) 2 - 4fc 2 (l + v + k 2 )} 



7T 



xcosech—P ) e" 



:(l+fc"). 



+4 / ;//.' 



a+ikC, 

00 ( fc2 - ^ 4 - 2?A; 3 tanhC) t 



;i + A; 2 ) 2 sinhffc 



(A.10) 



Now, before writing down the final form of ip^((, to) we evaluate the integrals 
with respect to k in the right hand side of the above equation. For this first we 
rearrange the integrand in Eq.(A.lO) by simple multiplication and call them 
J 3 , 7 4 , J 5 and I 6 as given below. 



- „ (k 2 - k 4 - 2ik 3 tanh () ikC , A „ . 

dk- — ; r~o7o ; m e tk \ (A.lla) 

(1 + k 2 ) 2 sinh \k ' V ; 

f°° ,, (1 - A; 2 - 2i/ctanhC) , 2 . . (i+fc 2 ) .. , / a111 n 

7 ' = L A (l + ^ si nh P + " + * (A ' llb) 

/*°° J7 (! - fc2 - 20; tanh C) iLl+^l a+ikC , A , , * 



2' 

oo 



The integral J3 can be evaluated by finding the residue of the integrand 

(fc (l+fc^inhT^ ° e ' fcC at the P° les fc = * of order two and at the sim P le P ole 
k = 2in. The results are given by 
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R es (k = i) = -(( - 2)sech(, (A.12a) 
i?es(A; = 2m) = ^= \J sechC,e~^ + ^= tanhC^/sec/j-Ce'^. (A.12b) 



While writing the first term in Eq.(A.12b), we have approximated the results 
obtained as a series in terms of suitable functions. Adding the above two 
residue values we obtain the value of the integral J 3 as 



■■2m 



-(C - 2)sech( + ^V-^Ce- " + ^tanhC 

(A.13) 



x V7 sech( e 



5C 
" 3 



In the case of the integral I4, the integrand possesses a second order pole at 
k — % and there is one more simple pole at k — (— i±y (^i+fp- -0 which 
is however out of the contour. Thus the residue for the integrand function 

+ „ + k 2 )e iH±pa + iK ig Qbtained ag 



i?es(A; = i) = -(2 + v( + 2va)sech(, (A.14) 
2 



and hence the right hand side of Eq. (A.14) represents the value of the integral 



The integrand of the integral I 5 namely (1 ~f 2 ;Sg h e iS2 ¥- La+iki admits an 
essential singularity at k — 0, in addition to a simple pole at k — ^^(-j ± 
(2n+i) 2 ~ ~) which is also out of the contour. The residue corresponding 
to the essential singularity is found by collecting the coefficient of | after 
expanding the above function. Thus we obtain 
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oo ( _ ir+ n' (f) 2n-l & n-l (& _ ira n' 

Res{k = 0)= 2^ * — ^ (C + oc) 

n,n'=0 



x(2 2 "- 1 - 1) + 
(C + «) 



(ra'!) 

'tt\2l2 ; 



-Br, 



(n!) 2 (n'!) 



2 ) 6 B n+ i (2 2„+i_ 1)[1 _ 2tanhC 



x 



X- 



(n + l)! 2 (n 
(C + «) 2 



(n' + l) (ra'+l)(n' + 2) 
(n + 1)(2 2 "+ 3 - 1 



+ 



(f) 4 & 4 £n + 2 

> + 2)!(ra + 3)! 



(n' + 2)! 

(C + «) 2 



l-2tanh C ^ 
(n + 3) 



(n' + 3)(ra' + 4) 



+ 



(A.15) 



where -B^s are Bernoulli numbers and b = ^y^. Thus the right hand side of 
Eq. (A.15) gives the value of the integral 1$. In the case of the integral Iq, 

the integrand ^r^^^^^- el +k )a+lkC - possesses a second order pole at k — i 
and an essential singularity at k — 0. In addition, we have a simple pole at 
k = (—2 ± \J (2n+\) 2 ~ ^) w hich is again out of the contour. The residues 
at the pole fc = i and at the essential singularity fc = are respectively found 
to be 



Res(k = i) 



-1 



TCV " 



1 — 2i>£ — 4f a)sec/iC, 



(A.16a) 



oo m 



Res(k = 0)= E- 

m,n,rt'=0j=0 



-l) n+n '(6 - l) n B„ +m (f) 2ri+2m & n+2m a n ' 
(n!)(n + 2m)!(n')!(n' + 2j)! 



x(C + a) 



n'+2j 



(n + 1 + 2m) (n + 2 + 2m) 



2(f) 2 & 2 (C + «) tanhC 



(n + 1 + 2m)(n + 2 + 2m) (ri + 1 + 2j) 



(A.16b) 



where E' n s are Euler numbers. We evaluated the value of the residue given 
in (A. 16a) using Mathematica. The value of the integral 1$ is the sum of the 
right hand sides of Eqs. (A. 16a) and (A. 16b). 



Now the value of ^ (1) (C*o) is found by combining Eqs. (A. 13), (A. 14), (A.15), 
(A. 16a) and (A. 16b) and the final form is written as 



* (1) (C,*o) 



80 



27V2 

7T 



64 



6v 2 



sech( e 2 + p=tanh (ysech( e ^ 
2v(l + v)( + v 2 + Aav - lj sech(. 



(A.17) 



While writing the above, we have dropped few higher order terms due to 
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smallness in values that appeared in the residue of the essential singularity. 



B Evaluation of the integrals in Eq. (75) using residue theorem 



In this appendix, we evaluate the integrals found in Eq. (75) using standard 
residue theorem as done in the previous two cases. Eq.(75) is written as 



* (1) (C,*o) 



I POO 



dk 



/oo 
d((l - k 2 
-oo 



it J -oo (1 + A; 2 ) 3 

7T 2 

+2ik tanh () [sin ( + (cosC — — ) tanh Qsech( 

8 

x{e l -^ a e iK -e- tk<: } + (l + v)sech( / d( ( 2 

\_J — oo 



7T 



x [sin ( + cos (]sech ( — — I d( ( z sech A (tanh( 

8 J-oo 



(B.l) 



It may be verified that the last integral J^d^ ( 2 sech 2 ( tanh£ in the right 
hand side of Eq.(B.l) on evaluation vanishes. For evaluating some of the inte- 
grals found in Eq. (B.l) we use the values of the integrals given in Eqs. (A. 6), 
(A. 8a) and (A. 8b) and also the value of the following integral. 



/« = 



C 2 tanhCe^C- 



(B.2) 



The integrand in Eq.(B.2) is found to be analytic everywhere except at the 
pole C = i(2n + l)|, n=0,l,2,... The residue of the function ( 2 tanh(e^ is then 



found to be — a 

8sinh(-|) 

integral I 8 is written as 



1 + 2 



sinh(f) 



+ 



2sinh a (f ) 



and hence the value of the 



/oo 77T 

Now, using the value of the integral given in Eq. (B.3), we evaluate the fol- 
lowing integrals found in Eq. (B.l) by integrating them by parts successively. 



1 + 2 



+ 



e2 (1 — e" 



-2n\ 



sinh(f) 2sinh 3 (f) v 



(B.3) 
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d(( 2 sech 2 (e ±l<: = 



71 



sinh | 



±(lT - )+ (i±f)fii 

V ^ 4 ; sinhf 



± 



vre ±i(l_ e T2^ 



sinh 



3 7T 



(B.4a) 



dCC 2 ^ec/i 2 CtanhCe ±lC : 



ITT 



sinh | 



2 sinh | 



i7r 2 e ± f(l -e^ 27r ) 
16 sinh 3 f 



(B.4b) 



On substituting Eqs. (B.4a) and (B.4b) in Eq. (B.l) we obtain 



* (1) (C,*c 



(1 + v)iir 3 
16 sinh | 



COsh^ 7T 

J — cosh — 



sinh 



2 \sinhf 



+ 1 



sech( 



x (1 + A; 2 ) 2 } (sec/i|(/3 + 1) + sec/i|(/3 - 1)) - ^6(1 - 6) 
/-, ,2n2 , * /*°° ,, (A; - A; 3 - 2i/ctanhC) 

x(i + ^ 2 ) 2 -Vj - 2 7-00*" (T+pj 3 " 

xe !fcC (sec/i|(l - fc) + sech^(l + k)j . (B.5) 



Before evaluating the integrals in Eq.(B.5) we rewrite the same appropriately. 
We then evaluate the integrals in Eq. (B.5) one by one by finding the values 
of the residues at the pole (k = i) at different orders. The residue for the 

integrand function (sec/i§(/3 + 1) + sec/if {(5 - 1)) Ml^_2|^M) e ^^+ifcC 
at the pole k = i of order three is found to be 



Res{k 



71 



- — [2(26 - 1)C + 4(26 - l)a + l)]sech(. 
16 



(B.6) 



Next, we find the residue for the function e 
simple poles k — % and k 



lil±^+ifct (l-fc 2 -2»fctanhp t fh 
fc(l+fc 2 )cosh/3§ dX Ule 



(2n+l) , 
(1-v) 



-i± 



l-v- 1 



>ll+1 \2 — 1) ( which is out of contour). 
Further, at k — there is an essential singularity, the residue of which is not 
shown here due its unwieldy form. Thus the residue for the above function at 
the first order pole k = % is given by 



Resik = i) 



tt(1 -26) 2 



[2(26 - 1)C + 4(26 - l)a - l]sech(. 



(B.7) 
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The residues for the function (sech^(k - 1) + sech\{k + 1)) ^^ff^f^e^ 
at the pole k — i of order three and at the simple poles k = i(2n + 1) + 1 and 
k = i(2n + 1) — 1 are written as 



Resik = i) = ^(2C- l)sech(, (B.8a) 
16 

x{[(2n + l) 2 -2i(2n + l)] 

-2itanhC[l + i(2n + 1)]}. (B.8b) 

1 V ; J 7T ^ [(2n + l) 2 -2 + 2*(2n + l)] 3 

x{[(2n + l) 2 + 2i(2n+l)] 
-2itanhC[-l + i(2n + l)]}. (B.8c) 

It can be verified that the residue of the function (sech*(P + 1) + sechKfi - 1) 

^ 1 ~ fc fc (" 1 2 ^fc2) nh ^ 6 ( k Hht " Sit k — i vanishes. On summing up the values of the 
residues given in Eqs. (B.6), (B.7) and (B8) (dropping higher order terms in 
Eqs. (B.8b) and (B.8c) while adding) we obtain the following expression for 
^"(C,to). 



b{2b - 1)C + {b(2a + b - 1) - a} + 



122825 

x{(cosh2C-sinh2C) [25 cosh 2C(523 cos C - 1512 sin C) + 289 
(19cosC + 8sinC) + 15(766co< -2393 sin C) sinh2C]}] . (B.9) 

As residues due to the singularities in the plane except along the imaginary 
axis, lead to secular terms in the solutions, we take into account only singu- 
larities along the imaginary axis that is at k — i, dropping residues due to the 
other two singularities at k — i(2n + 1) + 1 and k = i(2n + 1) — 1. Thus the 
first order correction V^CC t ) given in Eq.(B.9) finally becomes 



* (1)(C ' to) * 4(1 _26)2 i b ( 2b ~ !)C + { & (2« + b-l)-a}\ sech(. (B.10) 
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